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A variational method for many electron system is applied 
to momentum distribution calculations. The method uses a 
generating two-electron geminal and the amplitudes of the 
occupancies of one particle natural orbitals as variational pa- 
rameters. It introduces correlation effects beyond the free 
fermion nodal structure. 



The characteristics of condensed matter systems are 
due to the motion and correlation of the electrons |Q. 
The electron motion can be observed by Compton scat- 
tering with photons or by positron annihilation. Recent 
experiments |^^| indicate that the momentum density 
even in simple metals cannot be well represented by a 
single Slater determinant state. Instead, the momentum 
density has to be constructed from a correlated state with 
average occupancies m of single particle states in between 
and 1 0. In other words, at the independent particle 
level there are only N occupancies different from zero, 
and these are equal to one. When correlation is intro- 
duced, we have an infinite number of occupancies differ- 
ent from zero, even though most of them will presumably 
be very close to zero. The sum rule 



N 



(1) 



remains always true (the factor 2 is due to the spin). The 
purpose of the present work is to give a simple and effi- 
cient calculation method to estimate the occupancies rij. 
To simplify the problem, we will consider approximations 
which neglect the spin. Therefore the present discussion 
applies to non-magnetic systems. 

If the many body state is given by the wave function 
VP, the first-order density matrix p is defined by 



P (r,v') = N J de**(r,0*(r',0 



(2) 



The eigenfunctions ipi of p, introduced by Lowdin as nat- 
ural orbitals S arc the most suitable set of one-particle 
functions to use in this discussion 



p = 2} j n i \ip i >< ipi 



(3) 



The natural orbitals form an orthonormal basis set. An- 
other set of orbitals that are naturally associated with 
many-body functions are the generalized overlap ampli- 
tudes p]. These orbitals are, however, linearly dependent 
and a canonical orthonormalization of them yields the 
natural orbitals. 

The range of the the first-order density matrix p(r, r') 
in real space is a fundamental property of quantum me- 
chanical systems since it determines the degree of locality 
of the bonding properties [Q. 

The two-particle reduced density matrix a contains all 
the information to discuss two-particle interactions V 2 
HU and the total energy can be expressed as 



E[a] = (N/2) Tr(Ko-) , 



(4) 



K(r 1 ,r 2 )=H 1 (v 1 )+H 1 (r 2 ) + (N-l)V 2 (r 1 ,r 2 ) , (5) 

where Hi is the one-body part of the hamiltonian. If 'F 
is given by single Slater determinant, as in the Hartree- 
Fock approximation or in the density functional theory 
, then the energy is even determined by a one-particle 
density matrix p, such as p = p 2 (idempotency) . Re- 
cently, Goedecker and Umrigar (GU) Jl(| proposed to 
relax the p idempotency and to use a a natural orbital 
functional for a. The GU functional gives still a partic- 
ular importance to the individual electron picture. 

In the present work, an alternative method is explored. 
One considers the ansatz proposed by Blatt Jll|-|l3t , 

* = const 9l a + W*)a + Wi)] N / 2 | > . (6) 



where a + (ipi) are creation operators of an electron in the 
state ipi- Id coordinates space, \& is an Antisymmetrized 
Geminal Product (AGP) 



^ = const Det|0(r s ; — rj)\. 



(7) 



The generating geminal <f> has a diagonal expansion in 
the natural orbitals 



"i,r 2 ) = 3* ^*(ri)^(r 2 ) 



(8) 



In practice, the total energy becomes a functional 
E\gi,if)i]. Thus, gi and ipi are determined by minimiz- 
ing this functional. Such calculations have been done 
for some molecules |14 . The Stochastic Gradient Ap- 
proximation (SGA) optimization Jl5[ is particularly ap- 
propriate for the present problem since the variational 



1 



parameters can be determined avoiding the explicit de- 
termination of the total energy. 

The AGP is the N particle component of the the BCS 
state @. In the limit of N large, the AGP and BCS 
states become identical. If one set gi — Vi/v,i with 
\ui\ 2 + \vi\ 2 — 1, then m = \ui\ 2 . Therefore the elec- 
tron momentum distribution n(p) is given by the simple 
formula 11711 



i(p) = 2j2h\ 2 |<p| Vi> 



(9) 



The expectation value in the AGP of two-particle oper- 
ators can be found in ref. fjl3|| . 

For a 2 electron system, the present scheme is equiv- 
alent to a configuration interaction calculation and the 
two-particle reduced density matrix a is given by a pure 
state \<fi > 



a = <p >< 



(10) 



thus a = a 2 . 

The hydrogen molecule is a good example to illustrate 
the method. The bonding and antibonding orbitals are 

^ o(r) = [2(l + SW/2 [Mr) + Il{v)] ' 



Mr) = [2(1-S)]V2 [/fi(r) ~ /i(r)] ' (12) 



where f^i = \J a 3 /Tre~ a \ r ~ RR - L \ are Is atomic orbitals, 
S is the overlap integral (varying from to 1) and a is 
a variational parameter (varying from 1 to 1.66). The 
two-body wave function <j> can be approximated by 



(r,r') = g Mr)Mr') + 9iMr)Mr') ■ (13) 



Then 



So = -h=[1 



il/2 



V2 v/l + K 2 



1 
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fl- 



il/2 



V2 VTTk 2 ' 

where k is a function of S and of the integrals 



U=<f L f L I — | f L f L >, 



V =< f L f R | 1 f L fn >, 

ri2 



t =< fhfh I — I SlSr >, 



(14) 



(15) 



(16) 



(17) 



(18) 



J=<hh I — | f R f R > 

ri2 



(19) 



For large d, S w 0, k = 2t/(U - V"). Therefore, when 
d — * oo, go = —g\ = l/y/2 and 

<I>(t,i>) - g (f R (r)f L (r') + f L (r)f R (r% (20) 

This means that the correlation effects drive the elec- 
trons back on their own atoms like in the Heitler-London 
ansatz. 

For the linear chain molecule H4, the | ipi > (i = 
0, 1, 2, 3), have i nodes. In momentum space, < p | ipQ > 
is peaked at p = 0, but the < p | ijji > (for i = 1, 2, 3) are 
peaked at higher momenta. When d is small, only | ipo > 
and I tpi > are occupied, while, in the limit d — > 00, 
the SGA method yields go = g\ = —92 = —33- In the 
Hartree-Fock approximation [|18| , the momentum density 
n(p) of the chain H32 is more similar to that of a free- 
electron gas, with a given Fermi momentum pp , rather 
than that of the hydrogen atom. However, when the oc- 
cupation number can vary, one expects n(p) to develop 
high momentum tails. Recent experiments probing the 
electron momentum distribution in simple metals f^|J|] 
have observed similar tails. 

The Homogeneous Electron Gas (HEG) is another in- 
teresting limit for solids. In this system, the plane waves 
are the natural orbitals and the total energy per particle 
e = E/N is a function of the density parameter r s (i.e. 
the radius of the volume taken by one electron) . The dif- 
ference between the interacting and free HEG momen- 
tum densities for different r s yields the Lam-Platzman 
correction ]20| ] within the density functional theory. 

Csanyi and Arias plj computed the GU energy func- 
tional in the HEG and minimized the result with respect 
to the the occupancy n(k). At high density (small r s ), 
the result seems to reproduce the correct RPA limit and 
n(k) has a Daniel- Vosko like momentum dependence p2| . 
However when r s = 1, one finds e = 0.546 (a.u.), while 
the Diffusion Quantum Monte Carlo gives e = 0.596 
(a.u.) JT9|. This is a quite surprising result, since a 
variational result should be always greater than the ex- 
act energy. The reason is that the two-particle reduced 
density matrix a has been varied over too large class of 
functions: the restriction to iV-representable a has not 
been imposed. In other words, one cannot find a many- 
body state yielding this a. The AGP is by definition 
iV-representable. Therefore, it provides a general varia- 
tional scheme for many-electron system. When the AGP 
is applied to the interacting (with Coulomb repulsion) 
HEG, one finds the independent particle occupation p3| . 
However, correlation effects in the particle occupation 
may appear in an inhomogeneous electron gas. 

A crucial question is whether or not the fixed node ap- 
proximation, used by the Quantum Monte Carlo (QMC) 
simulations, gives a significant momentum density error 
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in realistic extended systems. A recent QMC calcula- 
tion for solid Li |^4j corrects about 30 % the discrepancy 
between the experimental Compton profiles and the 
density functional result. The fixed node approximation 
might be a cause of the remaining discrepancy. The Lam- 
Platzman correction, which is in the same nodal struc- 
ture, gives a similar result (see Fig. [I]). 




p(a.u.) 



FIG. 1. Total valence-electron Compton profiles of Li 
along (1 0). The solid line is the LDA calculation, the 
dotted line is the LDA with Lam-Platzman corrections and 
the dashed line is the AGP with |A(k)| = 0.1 a.u.. 

It is therefore worthwhile to investigate schemes be- 
yond the free fermion nodal structure like the AGP. In 
solid Li one can approximate the natural orbitals by the 
Kohn-Sham orbitals Q and do the following BCS ansatz 
for the occupation amplitudes [j25|,pG| 



.9(k) 



A(k) 



e(k) + ^e 2 (k) + A 2 (k) 



(21) 



The band energy e(k) is zero at the Fermi level and A(k) 
can be either calculated variationally or fitted to the ex- 
periment. Fig. ^ shows that important correlation effects 
can be observed in the Li Compton profile if |A(k)| is 
about 0.1 a.u.. 

In conclusion, the present paper presents a total energy 
functional of natural orbitals. The method goes beyond 
the Slater determinant nodal structure. For 2 electron 
systems, it is equivalent to a configuration interaction cal- 
culation. It can capture important correlation effects in 
the electron momentum density calculation. The knowl- 
edge of these effects is crucial for a proper interpretation 
of the experimental spectra. 
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